Non-Gaussian equilibrium in a long-range Hamiltonian system 



Vito Latora", Andrea Rapisarda'' and Constantino Tsallis^ 
°' Dipartimento di Fisica e Astronomia, Universitd di Catania, and INFN sezione di Catania, 

Corso Italia 57 / 95129 Catania, Italy 
* Laboratoire de Physique Theorique et Modeles Statistiques, Universite Paris-Sud, 91405 Orsay Cedex, France 
Centra Brasileiro de Pesquisas Fisicas, Rua Xavier Sigaud 150, 22290-180 Rio de Janeiro, Brazil 

(February 1, 2008) 

We study the dynamics of a system of N classical spins with infinite-range interaction. We show 
that, if the thermodynamic limit is taken before the infinite-time limit, the system does not relax 
to the Boltzmann-Gibbs equilibrium, but exhibits different equilibrium properties, characterized by 
stable non-Gaussian velocity distributions. Levy walks and dynamical correlation in phase-space. 
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' Though not always clearly stated, standard equilibrium thermodynamics [|l]-D is valid only for sufficiently short- 
' range interactions. This is not the case, for example, for gravitational or unscreened Coulombian fields, or for systems 
^ ' with long-range microscopic memory and fractal structures in phase space. The increasing experimental evidence 
'"^ . of dynamics and thermodynamics anomalies in turbulent plasmas Band fluids , astrophysical systems [p[-p^ , 

nuclei |13 1^ and atomic clusters granular media glasses 17, ij] and complex systems p9| , p0| found in the 
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last years, provide further motivation for a generalization of thermodynamics. 
^Kf-^ I In this paper we consider a simple model of classical spins with infinite range interactions pl|"p4[|, and we show that, 
^ . if the thermodynamic limit is performed before the infinite time limit, the system does not relax to the Boltzmann- 
' Gibbs (BG) equilibrium, but exhibits different equilibrium properties characterized by non-Gaussian velocity distribu- 
\l tions. Levy walks and dynamical correlation in phase-space, and the validity of the zeroth principle of thermodynamics. 
Our results show some consistency with the predictions of a generalized non-extensive thermodynamics recently pro- 
. posed 1^ , ^ . The Hamiltonian Mean Field (HMF) model describes a system of N planar classical spins interacting 
^-H ' through an infinite-range potential ||2l]]. The Hamiltonian can be written as: 
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'■^ where 9i is the ith angle and pi the conjugate variable representing the angular momentum (or the rotational velocity 
^ since unit mass is assumed). The interaction is the same as in the ferromagnetic X-Y model ||^, though the summation 
is extended to all couples of spins and not restricted to first neighbors. Following tradition, the coupling constant in 
the potential is divided by N. This makes H only formally extensive {V ^ N when N oo) p5|-p8|, since the energy 
1^ remains non-additive, i.e. the system cannot be trivially divided in two independent sub-systems. The canonical 
, analytical solution of the model predicts a second-order phase transition from a low-energy ferromagnetic phase with 
' magnetization M ~ 1 (M is the modulus of M = i where m.i — [cos{9i), sin{9i)], to a high-energy one 

5^ : where the spins are homogeneously oriented on the unit circle and Af ~ 0. The caloric curve, i.e. the dependence 
of the energy density U = E/N on the temperature T, is given hy U = + ^(1 — M^) and shown in Fig. 1(a). The 
critical point is at energy density Uc = 0.75 corresponding to a critical temperature = 0.5 pl| . The dynamical 
behavior of HMF can be investigated in the microcanonical ensemble by starting the system with water bag initial 
conditions (WBIC), i.e. 9i = for all i (M = 1) and velocities uniformly distributed, and integrating numerically the 
equations of motion [22| . As shown in Fig. 1(a), microcanonical simulations are in general in good agreement with 
the canonical ensemble, except for a region below Uc, where it has also been found a dynamics characterized by Levy 
walks, anomalous diffusion and a negative specific heat [Q. Ensemble inequivalence and negative specific heat 
have also been found in self-gravitating systems [ |j, nuclei and atomic clusters Jl3|-p^ , though in the present model 



such anomalies emerge as dynamical features |^,^. In order to understand better this disagreement we focus on 
a particular energy value, namely U = 0.69, and we follow the time evolution of temperature, magnetization, and 
velocity distributions. 

In Fig. 1(b) we report the time evolution of 2 < K > /N, a quantity that, evaluated at equilibrium, is expected to 
coincide with the temperature (< • > denotes time averages). The system is started with WBIC and rapidly reaches a 
metastable or quasi-stationary state (QSS) which does not coincide with the canonical prediction. In fact, after a short 
transient time, 2 < K > /N shows a plateau corresponding to a N-dependent temperature Tqss{N) (and Mqss ~ 0) 
lower than the canonical temperature. This metastable state needs a long time to relax to the canonical equilibrium 
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state with temperature Tcan = 0.476 and magnetization Mean = 0.307. The duration of the plateau increases with the 
size of the system: in particular we have checked that the lifetime of QSS has a linear dependence on N, see Fig. 1(c). 
Therefore the two limits t oo and iV — > cx) do not commute and if the thermodynamic limit is performed before the 
infinite time limit, the system does not relax to the BG equilibrium. This has been conjectured to be an ubiquitous 
feature in non-extensive systems [|5), but it has also been found for spin glasses When N increases Tqss{N) 
tends to Too = 0.380, a value obtained analytically as the metastable prolongation (at energies below Uc = 0.75) of 
the high-energy solution (M = 0). We have also found that [Tqss{N) - T^] oc N'^/^ and Mqss « iV^i/^, see Fig. 
1(d). At the same time we have checked that increasing the size, the largest Lyapunov exponent for the QSS tends 
to zero. In this sense mixing is negligible and one expects anomalies in the relaxation process |3^. The fact that 
Tqss converges to a nonzero value of temperature for N ^ oo means that, when N is macroscopically large, systems 
can share the same temperature, though this equilibrium is not the familiar one. All this amounts to say that the 
zeroth principle of thermodynamics is stronger than what one might think through BG statistical mechanics, since 
it is true even when the system is not at the usual BG equilibrium. We have checked the robustness of the above 
results by changing the level of accuracy of the numerical integration and by adding small perturbations. We also 
verified that the QSS has a finite basin of attraction, by adopting different initial conditions, as for example double 
water bag (DWBIC). In Fig. 2 we focus on the velocity probability distribution functions (pdfs). The initial velocity 
pdfs (WBIC or DWBIC) , reported in Fig. 2(a) , quickly acquire and maintain during the entire duration of the 
metastable state a non-Gaussian shape , see Figs. 2(b) and 2(c). The velocity pdf of the QSS is wider than a Gaussian 
for small velocities, but shows a faster decrease for p > 1.2. The enhancement for velocities around p ~ 1 is consistent 
with the anomalous diffusion and the Levy walks (with average velocity p ^ 1) observed in the QSS regime The 
following rapid decrease for p > 1.2 is due to conservation of total energy . The stability of the QSS velocity pdf can 
be explained by the fact that, for N oo, Mqss — * and thus the force on the spins tends to zero with N, being 
Fi = —MxSinOi + MyCosOi. Of course, for finite N, we have always a small random force, which makes the system 
eventually evolve into the usual Maxwell-Boltzmann distribution after some time. We show this for small systems 
(N=500,1000) at time t=500000 in Fig. 2(e). When this happens. Levy walks disappear and anomalous diffusion 
leaves place to Brownian diffusion . A possible frame to reproduce the non-Gaussian pdf in Fig. 2 (b) could be the 
non-extensive statistical mechanics recently proposed with the entropic index q ^ 1. This formalism provides, 

for the canonical ensemble, a q-dependent power-law distribution in the variables pi , 9i . This distribution has to be 
integrated over all 0i and all but one pi in order to obtain the one- momentum pdf, Pq{p), to be compared with the 
numerical one, Pnumip), obtained by considering, within the present molecular dynamical frame, increasingly large 
N-sized subsystems of an increasingly large M-system. Within the M >> N >> 1 numerical limit, we expect to 
go from the microcanonical ensemble to the canonical one (the cut-off is then expected to gradually disappear as 
indeed occurs in the usual short-range Hamiltonians), thus justifying the comparison between Pq{p) and Pnum{p)- 
The enormous complexity of this procedure made us to turn instead onto a naive, but tractable, comparison, namely 
that of our present numerical results with the following one-free-particle pdf ||2^ Pip) = [1 — iw^^^ ~ '?)p^]^^*-^~'^'' 
, which recovers the Maxwell-Boltzmann distribution for q = I. This formula has been recently used to describe 
successfully turbulent Couette- Taylor flow ||] and non-Gaussian pdfs related to anomalous diffusion of Hydra cells in 
cellular aggregates In our case, the best fit is obtained by a curve with q — 7 , T = 0.38 as shown in Figs. 2 
(b) and 2 (c). The agreement between numerical results and theoretical curve improves with the size of the system. 
A finite-size scaling confirming the validity of the fit is reported in panel (d), where A = Pth — Pnum, the difference 
between the numerical results and the theoretical curve for g = 7, is shown to go to zero as a power of N (for four 
values of p). Since g > 3, the theoretical curve does not have a finite integral and therefore it needs to be truncated 
with a sharp cut-off to make the total probability equal to one. It is however clear that, the fitting value g = 7 is 
only an effective non-extensive entropic index. Similar non-Gaussian pdfs have also been found in turbulence and 
granular matter experiments Pjl^, though this is the first evidence in a Hamiltonian system. In Fig. 3 we verify, 
through the calculation of the fractal dimension D2 , that a dynamical correlation emerges in the /i-space before 
the final arrival to a quasi-uniform distribution. During intermediate times some filamentary structures appear, a 
similar feature has recently been found also in sclf-gravitating systems |l^], which might be closely related to the 
plateaux observed in Fig. 1(b). We learn from the curves in Fig. 3(c) that, since they do not sensibly depend on N, 
the possible connection does not concern the entire /i-space, but perhaps only the small sticky regions between the 
"chaotic sea" and the quasi-orbits p^ . 

Metastable states are ubiquitous in nature. Their full understanding is, however, far from trivial. They basically 
correspond to local, instead of global, minima of the relevant thermodynamic energy. The two types of minima are 
separated by activation barriers which, at the thermodynamic limit, can be low, high or infinite, all of them presumably 
occurring in nature. The last case yields of course to quite drastic consequences. Moreover, the local minimum can 
either make the system to live in a smooth part of the a priori accessible phase space, or it can force it to live in a 
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geometrically more complex (e.g., multifractal) part of the phase space. The richness of such situation is what makes 
interesting the study of glasses, nuclei, atomic clusters, self-gravitating and other complex systems. It is natural to 
expect for such systems that the infinite size and infinite time limits are not interchangeable. What has emerged 
quite clearly here is that thermodynamically large systems with long-range interactions belong to this very rich class. 
We have verified that the usual attributes of thermal equilibrium: zeroth principle at finite temperatures, robustness 
associated with a finite basin of attraction in the space of the initial conditions, stable distribution of velocities, are 
satisfied, but they systematically differ from what BG statistical mechanics has make familiar to us along the last 
130 years. Our findings indicate some consistency with the predictions of non-extensive statistical mechanics [ p5[ , 
though a firm and unambiguous connection remains a challenge for future studies. In particular we believe all these 
features not to be exclusive of the present HMF model. Similar scenarios are expected for systems with say two-body 
interactions decaying like r~" for < a < etc , where ac is equal, for classical systems, to the space dimension p7| , p8| . 
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FIG. 1. (a) Caloric curve: microcanonical numerical results for N = 10000, 100000 are compared with equilibrium theory 
in the canonical ensemble. The dashed vertical line indicates the critical energy. Water bag initial conditions (WBIC) are used 
in the numerical simulations. Temperature is computed from T = 2 < K > /N , where < ■ > denotes time averages after a 
short transient time to = lO'^ (not reported here), (b) Microcanonical time evolution of 2 < A" > /N, for the energy density 
U = 0.69 and different sizes. Each curve is an average over typically 100 — 1000 events. The dot-dashed line represents the 
canonical temperature Tcan ~ 0.476. The quantity 2 < K > /N which starts from an initial value 1.38 (V^ = and K = U ■ N 
in WBIC), does not relax immediately to the canonical temperature. The system lives in a quasi-stationary state (QSS) with 
a plateau temperature Tqss{N) smaller than the expected value 0.476. Lifetime of QSS increases with N and the value of their 
temperature converges, as A'^ increases, to the temperature Too = 0.38, reported as a dashed line. Log-log plot of QSS lifetime 
(c) and Tqss{N) — Too (d) are reported as a function of the size N. Lifetime diverges linearly with N, and Tqss{N) converges 
to Too = 0.38 as iV"^/^ (see fit shown as a dashed line). Note that from the caloric curve one gets = r+l-2;7 = r-0.38. 
Therefore from the behaviour reported in panel (d), being Too = 0.38, one gets Mqss = N~^^^ . Results are similar when we 
consider double water bag initial conditions (DWBIC), i.e. 6i = G for all i and velocities uniformly distributed in (— P2, — Pi) 
and (pi,p2). In figure we report the case p\ = 0.8, p2 ~ 1.51. 
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FIG. 2. Time evolution of the velocity probability distribution function (pdf) for U = 0.69 and different sizes of the system, 
(a) At time t = we start with a single (WBIC) or a double (DWBIC) water bag velocity pdf. (b) In the transient regime 
where 2 < K > /N shows a plateau corresponding to Tqss{N) and the system lives in a Quasi-Stationary State (QSS), the 
velocity pdfs do not change in time and are very different from the Gaussian canonical equilibrium distribution (full curve). 
The pdfs at time t — 1200 for = 1000, 10000, 100000 show a convergence towards a non-Gaussian distribution which can be 
fitted by means of a power-law analytical curve (dashed curve) consistent with the generalized non-extensive thermodynamics 
[ psi proposed by Tsallis and characterized by q = 7 and T = 0.38, see text. The theoretical curve has been truncated with a 
sharp cut-off in order to have total probability equal to one, see text, (c) The same curves shown in (b) ai t — 1200 are reported 
in linear scale, (d) We show the difference. A, between the theoretical curve and the numerical results, as function of N for the 
four values of p indicated by arrows in panel (c). (e) We show the numerical pfds at t — 500000 for N — 500 and 1000. We get 
an excellent agreement with the Gaussian canonical equilibrium distribution at temperature T = 0.476. 
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FIG. 3. Correlation in /i-space. (a) We show the /x-space, i.e. the angle and momenta of the A'^ particles, for U — 0.69 and 
A'' = 10000 at different timescales, starting from WBIC. Although the initial configuration is uniform, structures emerge and 
persist for a very long time before dissolving again at equilibrium. A way to measure these correlations is by means of the 
correlation integral g^l C'(^) — ~ '^ij) where dij is the Euclidean distance between two points of the ^-space. In 

general C(r) — r^^, where D2 is the correlation fractal dimension, (b) By reporting the logarithm of C{r) vs the logarithm of 
r, a linear behavior over several decades is found. The fractal dimension thus extracted is reported in (c) vs time (an average 
over 50 events is considered). In the same time scale where we find the QSS, the correlation dimension is in between 1 and 2. 
The particles are fully spread in the /x-space only at equilibrium. As time increases, D2 grows continuously from 1 to 2. 
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